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Abstract. In this work, we consider a general class of nonlinear Volterra 
integro-differential equations with Atangana-Baleanu derivative. We 
use the operational matrices based on the shifted Legendre polynomials 
to obtain numerical solution of the considered equations. By approxi¬ 
mating the unknown function and its derivative in terms of the shifted 
Legendre polynomials and substituting these approximations into the 
original equation and using the collocation points, the original equation 
is reduced to a system of nonlinear algebraic equations. An error esti¬ 
mate of the numerical solution is proved. Finally, some examples are 
included to show the accuracy and validity of the proposed method. 


1. Introduction 

In recent years, several definitions for fixed and variable order fractional inte¬ 
gral and derivative operators have been introduced, such as Riemann-Liouville 
[7, 15, 30], Griinwald-Letnikov [30], Caputo [30], Caputo and Fabrizio [3], Atan- 
gana and Baleanu (AB) [1] (to see more, refer to [6, 15, 30, 27, 31]). 

Many problems in different fields of science and engineering such as viscoelas¬ 
ticity and damping, diffusion and wave propagation and chaos can be modeled by 
fractional derivative and integral operators [4, 5, 13, 15, 19, 26, 28, 29, 31]. For 
example, Pedro et al. [23] have investigated the drag force acting on a particle 
due to the oscillatory flow of a viscous fluid. 

Fractional integro-differential equations are used in different applications of 
physics, chemistry, biology, engineering and applied mathematics. For example, 
Zakes and Sniady in [32] have studied the dynamic behavior of multispan uniform 
continuous beam arbitrarily supported on its edges subjected to various types of 
moving noninertial loads. Li et al. have used differential equations to model a 
memory behavior of shape-memory polymer [16]. 

In most cases. It is not easy to obtain exact solution of most these equations. 
Therefore, some researchers have proposed several approximation and numerical 
methods for solving such equations. For example, Rakhshan and Effati applied a 
generalized Legendre-Gauss collocation method for solving nonlinear fractional 
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differential equations [24], Ganji et al. applied the fifth-kind Chebyshev poly¬ 
nomials to obtain solution variable orders differential equations [12]. Zhao and 
Wang solved space-time fractional partial differential equations using fast finite 
difference methods [33]. Liu et al. used a preconditioned fast quadratic spline 
collocation method for partial differential equations [17]. In [25], authors solved 
Fredholm integral equations by the reproducing kernel Hilbert space method. A 
numerical scheme based on the Bernstein polynomials for solving diffusion-wave 
equations has been proposed in [11] (to see more, refer to [7, 8, 14, 21]). 

In recent years, by increasing number of applications of fractional integro- 
differential equations, researchers attention have been given to use of orthogonal 
basis functions. The concept of operational matrices of these basis functions can 
be employed for the solution of the problems [20]. Recently, many researchers 
applied orthogonal basis functions for approximating. For example, Liu et al. 
solved variable order fractional differential-integral equations using Chebyshev 
polynomials [18]. Nemati et al. applied Legendre polynomials for solving two- 
dimensional nonlinear Volterra integral equations [22]. In [9] and [10], Jacobi 
polynomials for solving differential equations of variable order have been pro¬ 
posed. By using these basis functions and their operational matrices, the main 
problem reduces to an algebraic system, which can be solved. Then, the approx¬ 
imate solution can be calculated. 

In this work, we study the following type of the nonlinear Volterra integro- 
differential equations 

+ XTu{t), tG[0,l], . . 

1^(0) = iro, 

where 

Tu{t) = f K{t, s) (f>{s,u{s)) ds, 

Jo 

and 0 < w < 1, the parameter A is a real constant, AT, </> and / are given 
known functions, u{t) is an unknown function. denotes the Caputo AB- 

derivative operator. 

The organization of this paper is as follows: in Section 2, we present basic 
definitions of fractional calculus. Section 3 is introduced the main properties 
of the shifted Legendre polynomials (SLPs) and operational matrices based on 
the SLPs. Section 4 is devoted to proposing a new numerical method for solving 
problem (1.1). In Section 5, an error estimate is proved for the numerical solution. 
Section 6 includes some illustrative examples. Finally, we conclude the paper in 
Section 7. 


2. Basic definitions and theorems 
Definition 2.1 (See [30]). Let 0 < cu < 1. The RL-integral is defined as 


j {t- sY ^u{s) ds. 
^ Jo 
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The RL-integral of order cj satisfies the following property 


* r(r; + l + a;) 


Definition 2.2 (See [30]). Let 0 < w < 1, it G iL^(0,1) and AB{uj) be a nor¬ 
malization function suchthat AB[Qi) = AB{1) = 1 and AB{uj) = \ — uj + —— 


Then 


(1) The Riemann AB-derivative is given as 


ABRrB„,(+\ _ AB{uj) d 


^D^uit) = 


1 — uj dt 


where E^{t) = 


r(wi +1) 


E^{-- - {t - sY)u{s) ds, 

1 — CO 


is the Mittag-Leffler function. 


(2) The Caputo AB-derivative is defined as 


-'D'^u{t) = f E^{- ^ {t - sY)u{s) ds. 

^ ~ ^ Jo l — UJ 


(3) The AB-integral is given as 

A / n 1 — ^ \ ^ 


f {t — sY ^uY) ds. 
Jo 


Let Uto = . „ , ' and /3uj = . ^, —r, then we can rewrite (2.1) as 
AB[uj) AB[ijj)T[ijj) 

= a^u{t) -t- f3ujT{uj + lY^I^u{t). (2.2) 

It is easy to report that the AB-derivative and AB-integral satisfy the following 
properties 

(1) = 0, C G M, 


( 2 ) ^ 


(-l)^a;T(u + l) ^ 

(1 — tij)*r(a;i + u-|-1) ’ ’ 


(3) - —tY, 

1 — oj 1 — OJ 

(4) = C{a^ + P^tY, C G M, 

(5) {ai_j + I3cj{v + u! + l)B{v + l,uj + where B{-,-) is the 

Beta function, 


( 6 ) 


(7) = uY) - ^^(0). 

Theorem 2.1 (See [1]). Let (^[0,1] he the space of all continuous functions de¬ 
fined on [0,1] and f,g^ (^[O, Ij. Then the following inequalities can he established 

t^^DffY) - ^^^Dfg{t)\Y < 6\\f{t) - 5(t)||oo, 

t^^'Dffit) - ^^^DfgiYY < 5\\f{t) - 5(t)||oo, 
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where 5 is a constant number. 

Theorem 2.2. Suppose that f and g satisfy the assumptions of Theorem 2.1, 
then we have 

t^irf{t)-^^irg{t)\\^ < e\\f{t)-g{t)\U 
where e = + /3aj • 


Proof. According to definition of the AB-integral, we have 
t^Iffit) - (fit) - g{t)) lU 

= \\a^ (fit) - g{t)) + l3^T{uj + {f{t) - g{t)) ||oo 

< acoWfit) - git)\\oo + /3ujT{uj + (fit) - g{t)) ||oo 

< {auj + Pu) \\f{t) -5f(t)||oo- 

(2.3) 

Taking e = + /d^j, the proof is complete. □ 

Theorem 2.3. Let 0 < w < 1 and s' : [0,1] x M —)■ M 6e o continuous function 
such that there exits f satisfying 

\f{ 9 ,ui{t)) - f{g,U 2 {t))\ < C\ui{t) - U 2 it)\, (2.4) 

for all ui,U 2 G M. Then the initial value problem given by equation (1.1) has a 
unique solution on (^[0,1] 

(1) If {oco + Puj)Ci < 1? when X = 0, 

(2) If {oco + /3uj){Ci + IAIC 2 C 3 ) < I, when A / 0. 

Proof. Taking the AB-integral of both sides of equation (1.1) gets 
uit) = ^^Iffit, uit)) + X^^IfTu{t) + uo 

= aujfit, u{t)) + l3ujT{uj + l)^^Iff{t, u{t)) + Xa^Tu{t) (2.5) 
+ A/3^r(a; + l)^^IfTu{t) + uq. 

We define the operator Tl : (7(0,1] —> (7(0,1] by 

2flu{t) = a^fit, uit)) + /3^r(w + l)^^Iffit, uit)) + X/d^Tiu + l)^^IfTuit) 

+ XotijjTuif) P (7, 

where c = uq — a^ifiO, uq). By (2.5), finding a solution of equation (1.1) in (7[0,1] 
is equivalent to finding a fixed point of the operator Tt. For ui,M 2 £ (7(0,1] and 
all t £ [0,1], we have 

|Tlni(t) - 2fiu2it)\ = \a^ (/(t,ui(t)) - /(t,U2(t))) + Xa^jT (ui(f) - U2it)) 

-)-/3^r(w l)^^If ifit,uiit)) - fit,U2it))) 

+ Xfl^Tiu, + l)^^IfT (ui(t) - U2it)) I 

< aco\fit,Uiit)) - fit,U2it))\ |A|aa;|F (ui(t) - U2(t)) I 

+ flu^Tiuj + l)\^^If ifit,uiit)) - fit,U2it))) I 

+ |A|/?^r(a; + l)\^^IfT (ui(t) - U2(t)) |. 
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?3 = sup K{t,s), (2.6) 

t,se[o,i] 

and according to (2.4), we have 

\f{t,Ui{t)) - f{t,U2{t))\ < ?l|ui(t) - U2{t)\, 

\(l){t,Ui{t)) - (l){t,U2{t))\ < ^2\ui{t) - U2{t)\. 

In veiw of (2.6), (2.7) and according to Tu{t), we have 

|01ttl(t) - 0IU2(t)| < - U2it)\ + - U2{t)\ 

+ |A|aa;?2?3|'Wl(i) “ U2it)\ + |A|1Ul(t) - U2{t)\ 

= (aco + /3a;)(?l + |A|?2?3 )||w 1 “ ll2||oo. 

Therefore the operator Tl is a contraction. The statement follows from Banach’s 
Fixed Point Theorem. □ 

3. The SLPs and operational matrices based on the SLPs 

3.1. The SLPs. The SLPs on [0,1] consist the set of orthogonal functions. The 
analytic form of the SLP of degree n is defined by 

n 

£o(t) = l, Ci{t) = 2t-1, £„(^) = ^a„,fc^^ (3.1) 


where 


{n — k)\{k\y 

For two arbitrary functions f,g£ L^[0,1] the inner product and norm are defined, 
respectively, by 


(_l)n+fc(^^^)! 


if, 9 )= [ f{t)g{t)dt, 
Jo 


ll/lli(o,i) = vW)- 

Any arbitrary function u G L^[0, 1] can be expaned using the SLPs as follows 


i{t) = '^UiCi{t) = 


where 


Ui = Q ^ [ u{t) Ci{t) dt, 
Jo 


Q = [qij] bj = o,i,---,M, 

and the elements qij given by 


' 10 , i^j. 
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The function u{t) can be approximated by taking only the first M + 1 terms in 
(3.2) as 

M 

u{t) ~ UM{t) = ^ = u'^ £(i), (3.3) 

i=0 

where U = [uo,ui, and 

£(t) = [£o(t),A(t),-- - (3.4) 

Similarly, any function Kit, s) in T^([0,1] x [0,1]) can be expanded in terms of 
the SLPs as 


K(t,s) ~ £^(t)/C£(s), 
where K, is an (M + 1) x (M + 1) matrix as 

^ = [ki,j]’ b j = 0,1,... ,M, 
and the shifted Legendre coefficients kij are qiven by 

((iL(t,s),£(t)),£(s)) 


kij — 


I A: 


\Cj{s) 


i,j = 0,1,. 


,M. 


3.2. Operational matrices of the SLPs. 

(1) The integration of the vectors £{t) defined by (3.4) can be approximately 
given by 

£(s) ds ~ P£(t), (3.5) 

where P is an (M + 1) x (M + 1) matrix which is called the operational 
matrix of integration based on the SLPs and is given by 



P = - 
2 


1 1 0 

=1 0 ^ 

3 3 


0 0 0 
0 0 0 


0 

0 

-1 

2M-1 

0 


0 

-1 

2M+1 


0 

0 


2M-1 

0 


(2) It is needed to evaluate the product of £(t) and £^(t), that will give the 
product operational matrix of the shifted Legendre basis. First, we write 


i+j 

Ci{t) Cj{t) = ^ Or Cr{t)- 

r=0 


(3.6) 


The coefficients Ur, r = 0,1, ■ ■ ■ ,i+j, are computed in the following man¬ 
ner. Multiplying both sides of equation (3.6) by Cm{t), m = 0,1, • • • , i+j, 
and integrating the result yields 


/ Ci{t)Cj{t)Cm{t)dt = '^ar / Cr{t) Cm{t)dt = - -—( 

Jo ^=0 -^0 2m + i 


(3.7) 
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Using (3.7) and the analytic form of the SLPs given by (3.1), we obtain 


am = {2m + l)[ Ci{t) Cj{t) Cm{t) dt 
Jo 


I j m 


{2m + 1) EEE ^ 2 ,/c t dt — {2m 


k=0 1=0 S=0 


where 


I j m 


^ Ci^k Cj^l Cm,s 

L. fc + / + s + i- 


k=01=0 S=0 


By substituting am into equation (3.6), we have 

i+j 

Ci{t)Cj{t) = Y, {2m + l)Ajj^m. (3.8) 

m=0 

Suppose that 

C = [co, Cl, • • • , Cm] , 

then using (3.8), we can write 

£(t)i2^(t)C~C£(t), (3.9) 

where C is an (M + 1) x (M + 1) matrix which is called the product 
operational matrix of the SLPs basis vector, and is given by 


C = [ci,j], bJ = 0,1,--- ,M, 


where 


2j + l 




Also, for an (M + 1) x {M + 1) matrix D = [dij], i, j = 0,1, • • • , M, we 
have 

£^(t)D£(t) ~ D*£(t), (3.10) 

where D* = [d*], i = 0,1, • • • , M is an 1 x {M + 1) vector given by 

M M 


di — {2i + 1 ) ^ ^ ^ ^ ^m,n,i dm,r 


m=0 n=0 

(3) The operational matrix of Rl-integral of order oj for vector £(t) is as 

^^I‘^S{t) ~ T“£(t), (3.11) 

where 

^ — [^i,j]^ ^i,j — "Y/ Ppd^i>PY{p ^ 1) ’ b J = • • • ) 


Ppd — 


Pj + dEu 


{-ib‘0 + 0! 


u - i){i'-y{p+i+uj + 1 )' 


with 
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The matrix is called the operational matrix of RL-integral based on 
the SLPs [2] . 

(4) The operational matrix of AB-integral of order oj for vector £(t) is ob¬ 
tained as 

= a^£{t) + 

According to (3.11), we have 

a^t) + P.riu + l)T-£(t) ^ 

= («a;/+/3a;r(u; + 1)T"^) £(t) = 

where I is an (M +1) x (M -|-1) identity matrix and -|- /3i^T{uj + 

1)T‘^. The matrix J‘^ is called the operational matrix of AB-integral 
based on the SLPs. 


4. Numerical method 


In this section, we introduce a numerical method for the solution of nonlin¬ 
ear Volterra integro-differential equations of the form equation (1.1). For this 
suppose, first we approximate as 



~ f 7 ^£(t). 

(4.1) 

By taking the AB-integral of (4.1) and using initial condition, we have 



u{t) ~ + uq. 

(4.2) 

We approximate uq as 

no — hf" £^{t). 

(4.3) 

By putting (4.3) in (4.2), 

we can write (4.2) as 



u{t) ~ 7£(f), 

(4.4) 

where 7 = U'^ J‘^ + . 

Let g{t) = f{t,u{t)) and h{t) = (p{t,u{t)). 

Then we 

approximate g{t) as 

g{t) ~ G^£(f). 

(4.5) 

For approximating Tu{t), first K{t^ s) and h{t) are approximated as 



K{t, s) ~ Sy’"{t) JC T(s), 

(4.6) 


h{t) ~ 

(4.7) 


By using (4.6) and (4.7), Tu{t) is approximated as 

ru(t)~ f £F{t))CQ{s)SF{s)Hds 
Jo 

= £^it)IC [ 2{s)£^{s)Hds 

Jo 

= £^{t)IC [ H£is)ds 

Jo 

= £^{t)ICH [ 2{s)ds = £^{t)ICHP£{s). 

Jo 


(4.8) 
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Substituting (4.1), (4.5) and (4.8) in equation (1.1), we have 

C/^£(t) = G^^{t) + S.^{t)lCHP^{s). (4.9) 

Let A = 1C HP. Using (3.10) to define the vector A*, and using this approxima¬ 
tion in (4.9) we obtain the following system 

-G^ -A* = 0. 

In veiw of (4.4), we rewrite (4.5) and (4.7) as 

/(t,7£(t)) = G'^£(t), 

By putting M + 1 points tj, i = 1, 2,..., M -|- 1 in equation (4.11), gives 

/(ti,7£(L)) - = 0, 

= 0 , 


(4.10) 

(4.11) 


where 


i 

M + 2' 


Equations (4.10) and (4.12) form a system of 3(M -|- 1) nonlinear equations of 
the vectors of U, G and H. By solving this system, the unknown parameters of 
the vectors of t/, G and H are obtained. Finally the approximate solution can 
be computed by (4.4). 


5. Error estimation 

In this section, we give an estimate for the error of the numerical solution of 
equation (1.1) with initial condition obtained by the proposed method in Section 

4. 


Theorem 5.1. Suppose that u and um satisfy the assumptions of Theorems 2.1 
and 2.2. Assume that u G 1] is the exact solution of equation (1.1) and 

uuit) = 7 ^il(t) is its approximation given by the method proposed in Section f. 
Then, we have 


\u{t) - UM 


2 < 


{M+ 1)12“^^+^' 


\ABCrMJ 


Vfu{t) - 


oo ^ 


6p 


I ABjuj 


ffu{t) - 


\/mT1M!22A^+i ’ 

ep 


yjM 1M!22A^+1 ’ 


where p = supggjo^i] |rt*'^“'"^^(0)|. 


Proof. Let that u*{t) is the interpolating polynomials to u{t) at points ti, 
i = 0 , 1 ,..., M are the roots of the shifted Chebyshev polynomials of degree 
M + 1, then we can write 


u{t) — u*{t) 


rt(^+i)(6') 
(M + 1)! 


M 


]^(i ti), 

i=0 


where 6 G [0,1]. 
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According to Chebyshev interpolation nodes, for t G [0,1] it gives 


\u{t) — = 


P 


(M + l)!22Af+i’ 

where p = sup^igjg ;^] 

Since the best approximation of a given function u G 1] in the finite 

subspace Sm = {C-oit)■, Ci{t),... is unique, we can write 

lk(t) - UM{t)\\l < - U*{t)\\l 

= [ \u{t) — u*dt 


(5.1) 


P 


dt = 


Jo \{M + l)\2‘^^+\ 

Taking the squared root of both sides of (5.1) gets 


P 


(M + l)!22^+i 


\u{t) - UM{t)\\2 < 


p 


(M + l)!22A^+i’ 

We know that ||.||oo < ^/M + lIMb- Then, in veiw of Theorems 2.1 and 2.2, the 
proof is complete. □ 


6. Test examples 

In this section, we consider several examples to confirm the accuracy and ap¬ 
plicability of the proposed method. 

Example 6.1. Let f{t,u{t)) = 2t — 1 — 2e + 2e^“*+*^ + (Er/f[|] 

—Erfi[^ — t])), A = 1, m( 0) = 1 and Tu{t) is given as 

Tu{t) = I ds. 

Jo 

Erfi(-) is the imaginary error function. The analytical solution of this example 
is u{t) = l — t + when oj = 1. We solved this example by the proposed method. 
By considering M = 5, the numerical results are shown in Figure 1 by different 
values of w. As seen in Figure 1, when the value of a; tends to 1, the approximate 
solution converges to the analytical solution. 

Example 6.2. Let f{t,u{t)) = e*(l -t- e^) — (l + e^*(—1-|-2t)) — 

A = 1, u(0) = 1 and Tu{t) is given as 

Tu{t) = f t^svJ{s)ds. 

Jo 

The analytical solution of this example is u{t) = e* when oj = 1. This example is 
solved by the proposed method with different values of ui. By considering M = 5, 
the numerical results are plotted in Figure 2. 

Example 6.3. Let f{t,u{t)) = 1 — A = 0 and tt(0) = 0. The analytical 

- 1 

solution of this example is u(t) = —- when a; = 1. By applying the proposed 

c 1 

method, we solved this example. By considering M = 5, the numerical results are 
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Figure 1. (Example 6.1) The exact and approximate solutions 
given by different values of u and M = 5. 



Figure 2. (Example 6.2) The exact and approximate solutions 
given by different values of u and M = 5. 


plotted in Figure 3 by different values of w. As seen in Figure 3, when the value 
of CO tends to 1, the approximate solution converges to the analytical solution. 


Example 6.4. Consider the following system of nonlinear Volterra integro-differential 
equation 


ABCD<^u(^t) = fi{t,u{t),v{t)) + XiTi{u,v){t), 
ABCD‘^v{t) = f 2 {t,u{t),v{t)) + X 2 T 2 {u,v){t), 
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Figure 3. (Example 6.3) The exact and approximate solutions 
given by different values of u and M = 5. 


where 

Ai = —1, A 2 = 1, 

fi{t,u{t),v{t)) = 1 + + 

, f 2 {t,u{t),v{t)) = 1-^- 

Ti{u,v){t) = /o +v{s))ds, 

T 2 {u,v){t) = f^tu(s)v(s) ds, 

.«(0) = 0, i;(0) = l. 

The analytical solutions of this example are u{t) = t and v{t) = 1 + t when a; = 1. 
By using the proposed method, we solved this example. By considering M = 5, 
the numerical results are reported in Figure 4, Table 1 and 2 by different values 
of Lo. The exact and numerical solutions obtained by different values of u are 
plotted in Figure 4. Furthermore, in Table 1 and 2, the exact and numerical 
solutions at some grid points with different values of oj are seen. 




(a) (b) 

Figure 4. (Example 6.4) (a) The exact and approximate solu¬ 
tions {u{t)) (b) The exact and approximate solutions {v{t)) and 
M = 5. 
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Table 1. (Example 6.4) Numerical results u{t) by different values 
of uj and M = 5. 


t 

UJ = 0.85 

UJ = 0.9 

UJ = 0.99 

UJ = 0.999 

UJ = 1 

0.1 

0.30231 

0.23076 

0.11235 

0.10123 

0.10000 

0.3 

0.51480 

0.44093 

0.31363 

0.30136 

0.30000 

0.5 

0.69220 

0.62946 

0.51305 

0.50131 

0.50000 

0.7 

0.83586 

0.79789 

0.71085 

0.70109 

0.70000 

0.9 

0.93217 

0.93812 

0.90640 

0.90066 

0.90000 


Table 2. (Example 6.4) Numerical results v{t) by different values 
of oj and M = 5. 


t 

UJ = 0.85 

UJ = 0.9 

UJ = 0.99 

UJ = 0.999 

UJ = 1 

0.1 

1.30285 

1.23099 

1.11236 

1.10123 

1.10000 

0.3 

1.52586 

1.44615 

1.31384 

1.30138 

1.30000 

0.5 

1.74109 

1.65396 

1.51429 

1.50142 

1.50000 

0.7 

1.97079 

1.86860 

1.71497 

1.70148 

1.70000 

0.9 

2.22022 

2.09479 

1.91644 

1.90162 

1.90000 


7. Conclusion 

In this work, we presented a numerical method based on operational matrices to 
solve Volterra integro-differential equations. For this do, we used the operational 
matrixces of integral, produce, fractional integral of Riemann-Liouville type and 
obtained the operational matrix of fractional integral of Atangana and Baleanu 
type. Then, by substituting these matrices and collocation points, we converted 
the Volterra integro-differential equations to an algebraic system. By solving 
this system, we obtained the approximate solution. Finally, some examples are 
considered to confirm the accuracy of the proposed method. 
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